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Abstract The question whether the linear absorption 
spectra of metal clusters can be interpreted as density 
oscillations (collective "plasmons" ) or can only be under- 
stood as transitions between distinct molecular states is 
still a matter of debate for clusters with only a few elec- 
trons. We calculate the photoabsorption spectra of Na2 
and Na^ comparing two different methods: quantum 
fluid-dynamics and time-dependent density functional 
theory. The changes in the electronic structure associ- 
ated with particular excitations are visualized in "snap- 
shots" via transition densities. Our analysis shows that 
even for the smallest clusters, the observed excitations 
can be interpreted as intuitively understandable density 
oscillations. For Naji", the importance of self-interaction 
corrections to the adiabatic local density approximation 
is demonstrated. 
PACS: 36.40.Vz, 31.15.Ew 



1 Introduction 

Among the earliest experiments providing insight into 
the electronic structure of metal clusters were measure- 
ments of the linear photoabsorption spectra To- 
day, still, they are among the most powerful probes of 
a cluster's structure, and understanding the effects of 
collectivity in absorption spectra in general is also one 
prerequisite for understanding the nonlinear regime that 
is being probed experimentally with increasing sophisti- 
cation 1^. In particular for positively charged sodium 
clusters, linear spectra have been measured for a broad 
range of cluster sizes and temperatures [Q. Typically, 
the experiments show a few strong absorption lines that 
exhaust most of the oscillator strength. However, despite 
the fact that the absorption spectra have been known for 
a long time, their theoretical interpretation is still being 
discussed. Since sodium is the nearly-free electron metal 
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par excellence, allowing to study metallic behavior in 
one of its purest forms, a lot of theoretical work over the 
years has been devoted to photoabsorption in Na clusters 
|[|,|,|,|[|l^,|ll|,[2[[l3[|l|,|l|. From these studies, two 
different and somewhat opposing points of view on the 
interpretation of the observed resonances emerged. On 
the one hand, small clusters can accurately be described 
in the language of quantum chemistry, understanding 
the excitations as transitions between distinct molecular 
(electronic) states. On the other hand, the experiments 
found an early, intuitive interpretation in terms of col- 
lective oscillations of the valence electron density against 
the inert ionic background, similar to the plasmon in 
bulk metals or the giant resonances in atomic nuclei. 
Since the strong delocalization of the valence electrons 
that characterizes nearly free electron metals is found 
even in the smallest Na clusters, it has been argued that 
the second interpretation should also be applicable to 
small clusters. 

A theoretically well founded and practically tested 
j5[^,|l^,|l2, 13; 1^ method for the theoretical investiga- 
tion of excitations in metal clusters is time dependent 
density functional theory (TDDFT) at the level of the 
adiabatic, time-dependent local density approximation 
(TDLDA). A refinement of TDLDA corrects for the self- 
interaction error leading to the scheme of a time-dependent 
self-interaction correction simplified by a global averag- 
ing procedure (TDSIC) ||l2|. A somewhat simpler, yet 
powerful, alternative is quantum fluid-dynamics. In an 
extension of earlier works 1^,0, quantum fluid-dynamics 
in a local current approximation (LCA) was recently de- 
rived making direct use of the ground-state energy 
functional of density functional theory. In this work we 
are comparing these methods using two very small clus- 
ters as test cases with a threefold aim: First, compar- 
ing LCA to TDDFT for exactly the same system allows 
to judge on the reliability of LCA results. Second, by 
comparing TDLDA to its on-average self-interaction cor- 
rected counterpart, we check the impact of self-interaction 
corrections on low-energy photoabsorption spectra. Third, 
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the combination of methods aUows us to demonstrate 
that, indeed, the experimental spectra for even the small- 
est clusters can be interpreted as valence electron density 
oscillations, leading to an intuitive understanding of the 
experimentally observed effects. 

In section |^ our theoretical methods are reviewed. 
Section || presents the results for Na2 and Na^, which 
are discussed and summarized in section ^ 



2 Theory 

Starting point of our investigations is the usual ground- 
state energy functional 

E[n;{R}] = T,[{ip}] + E^,[n] + J n{r)Vion{r; {R}) d\ 
°2 r /•n(r)n(r') ,3 , ,3^ 



■ d-^r' d-^r 



r - r' 
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R-i ~ Rj I 



(1) 



for a cluster of N ions of valence Z (for Na, Z — 1), 
valence-electron density n and ionic coordinates {R}. 
The noninteracting kinetic energy Ts is calculated from 
the Kohn-Sham orbitals {</?} and -ExcH denotes the ex- 
change and correlation energy for which we use the LDA 
functional of Ref. Generalized gradient approxi- 

mations, e.g. in general lead to a better descrip- 
tion of correlation effects in small systems. However, 
in the present case LDA is not too bad an approxima- 
tion due to the strong delocalization of the valence elec- 
trons. Vion is the sum of pseudopotentials Vion(r; {R}) = 
X^iLi Wps(|r — R,;|). We employ the smooth-core pseu- 
dopotential of Ref. |0. In combination with LDA it 
provides accurate bond lengths, which are important for 
optical absorption spectra and polarizabilities jlj,^. 
The ionic coordinates were obtained by self-consistent 
minimization of the functional ([|) with respect to both 
n and {R} and are given in Ref. The Kohn-Sham 
equations are solved directly, i.e., without basis sets, on a 
real space grid. We have verified that the coordinates ob- 
tained in the symmetry restricted optimizations of Ref. 
do not change noticeably if the optimization is done 
fully three-dimensional. In order to be consistent, the 
ionic configurations were reoptimized on the SIC level 
for the TDSIC calculations, as discussed below. 

Eq. (||) is also the key ingredient for the quantum 
fluid-dynamical LCA. A detailed discussion of the the- 
ory can be found in Ref. Therefore, we here re- 
strict ourselves to a brief sketch. The essential idea of 
LCA is to describe excitations as harmonic density os- 
cillations. The oscillating density is obtained from the 
scaling transformation 



where a(t) oc cos ujt and Sn is the so called density scal- 
ing operator 



Sn = (Vu(r)) + u(r) • V, 



(3) 



which contains - hallmark of a fluid-dynamical descrip- 
tion ~ a displacement field u(r). A similar, consistent 
transformation is also applied to the Kohn-Sham orbitals 
from which the density is constructed. From a variational 
principle and Eq. (Q), a set of coupled, partial differential 
eigenvalue equations for the Cartesian components of u 
is derived. The eigenvalues are the excitation energies, 
and from the solutions Ui,, absolute oscillator strengths 
and intrinsic current densities 



i^{r,t) = a^{t) u^{r) n^{r, a^{t)) 



(4) 



n(r,a(0) =e~"(*)^"n(r). 



(2) 



associated with a particular (the v-th) excitation are 
derived. Eq. (^ is the reason for the name "local cur- 
rent approximation". If a mode jj^ can be excited by the 
dipole operator D = —ez we call it a z-mode (or x,y, 
respectively) . 

It is important to note that the LCA is not a (semi) 
classical but a quantum mechanical method in the sense 
that it is derived on the basis of the quantum mechanical 
Kohn-Sham energy functional, which contains informa- 
tion on the quantal single-particle states in the kinetic 
energy. But the range of validity for the LCA is hard 
to assess formally [^. However, gathering experience 
on its performance, as done for earlier versions which 
were truly semiclassical methods |^ or approaches using 
a well guessed expansion basis of local operators , will 
lead to a better understanding. It is therefore one aim 
of the present work to test the accuracy of the LCA by 
comparing it to the well established TDLDA. 

For TDLDA and TDSIC, the numerical solution of 
the Kohn-Sham equations is done on a spatial grid with 
Fourier transformation for the definition of the kinetic 
energy. Accelerated gradient iteration is used for the 
static part. The dynamic propagation is done with the 
time-splitting method. For details of the technology see 
the review . The spectra are computed as described in 
22, 2^, We start from the electronic and ionic ground 
state configuration. An instantaneous boost of the whole 
electron density initializes the dynamical evolution ac- 
cording to time-dependent LDA or SIC. The emerging 
dipole momentum as function of time is finally Fourier 
transformed into the frequency domain. This delivers the 
spectral distribution of dipole strength. The initial boost 
is kept small enough for the method to produce the spec- 
tra in the regime of linear response. 



3 Results 

With the methods described in the previous section we 
first investigated the sodium dimer. At first glance, it 
could be expected that the two-electron system Na2 is 
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Table 1 Dipole excitations up to 4 eV for Na2. Energies EE in eV and oscillator strengths OS as percentages of the dipole 
sumrule mi — e^fi^ N Z / (2m) for LCA (superscript a), TDLDA (superscript b) and TDSIC (superscript c). Columns labeled 
"mode" indicate the direction of oscillation, see text for discussion. For comparison, we also list TDLDA results (superscript 
d) and experimental values (superscript e) from Ref. — indicates that the corresponding mode is not found in LCA, - 
that the strength in TDSIC was beyond numerical accuracy, no v. that no corresponding value has been given in the literature. 



EE" 


OS" 


Mode" 


EE' 


OS' 


Mode' 


EE'^ 


OS"^ 


EE'* 


OS"* 


EE-^ 


1.93 


30.9 


z 


2.09 


31 


z 


2.13 


36 


2.09 


31.4 


1.82 


2.56 


58.9 


x/y 


2.63 


56 


x/y 


2.65 


57 


2.52 


53.1 


2.52 


3.93 


1.7 


z 


3.67 


1 


z 


3.89 




3.28 


<1 


3.64 








3.72 


3 


x/y 


3.95 




no V. 


no V. 


no V. 



not described accurately in the quantum fluid-dynamical 
LCA. However, as seen from Fig. 0, TDLDA and LCA 
give similar results. Since the LCA currents were calcu- 
lated with the LDA functional, we first compare them 
to TDLDA and discuss TDSIC results later. It is im- 
portant to note that our TDLDA and LCA calculations 
were performed on exactly the same basis, i.e, using the 
same internuclear distance and pseudopotential. For a 
closer inspection we give in Table |l| the excitation en- 
ergies and percentages on the dipole sumrule for the 
excitations up to 4 eV that carry most of the oscilla- 
tor strength. Comparing columns 1-3 to columns 4-6 
reveals some noteworthy differences between LCA and 
TDLDA. First, for the lowest excitation, LCA gives an 
energy lower than TDLDA with a difference larger than 
the numerical uncertainty. Since LCA rests on a varia- 
tional principle, the fact that it leads to a lower excita- 
tion energy than TDLDA points at that it can be seen 
as an independent method. In this context we also note 
that our TDLDA for the z-excitation is consistent with 
the result in Second, whereas LCA is very accurate 
for the two low-lying, strong transitions, it does not seem 
to perform as well for the higher lying excitations. For 
technical reasons, the oscillator strength for the weak 
transitions is hard to asses in our TDLDA and TDSIC 
and could only be estimated in TDLDA. However, Fig. 
|l| shows that in comparison to TDLDA, the strength of 
the third peak is underestimated in LCA, and the LCA 
eigenvalue is too high. Tab. |l| shows why the TDLDA 
spectrum in Fig. |l| looks better. Whereas LCA only leads 
to one z mode, TDLDA around 3.7 eV gives excitations 
in both z and x/y direction. A "double-peak" structure 
has also been found in other TDLDA calculations ijl^ . 
Thus, we conclude that LCA gives some of the strength 
carried by transitions at higher energies, but it does not 
provide the same resolution as TDLDA. This is under- 
standable since TDLDA embraces the whole fragmenta- 
tion into the various one-particle-hole (Ip/i) states of the 
excitation spectrum, while LCA is bound to a "collective 
deformation path" . 

For Na2, TDSIC leads to overall similar results as 
TDLDA, therefore we do not show the spectrum in a 
separate plot. Besides the fact that the TDSIC spectrum 
does not show the "cut" in the low energy shoulder of 




Fig. 1 Experimental |g5|, LCA and TDLDA photoabsorp- 
tion spectrum S of Na2 in arbitrary units against excitation 
energy in eV. The line broadening is chosen phenomenologi- 
cally to match the experiment. 



the excitation at 2.09 eV, the main difference is that 
TDSIC leads to slightly higher excitation energies than 
TDLDA (see Tab. 0). The deviation is less than 0.05 eV 
for the low-lying, strong transitions, but it is more than 
0.2 eV for the higher ones (placing them at a similar en- 
ergy as LCA). The reason is that the bonding distance is 
0.1 ao smaller than in LDA |l^. This slight compression 
leads to a small blue shift. The higher excitations are 
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more sensitive because they are dominantly Iph transi- 
tions, and it is known that SIC has stronger effects on 
the single-particle states. A recent, detailed discussion of 
how several other approximations for E^c influence the 
dimer spectrum can be found in p5[ |. 

An intuitive understanding of the observed excita- 
tions can be obtained by looking at "snapshots of the 
density change" . In LCA, these are easily accessible since 
the local currents, Eq. m), obey the continuity equation 



at 



(5) 



Thus, one only needs to numerically calculate and then 
plot the divergence of ji, to obtain a visualization of the 
density change associated with the ly-th excitation at one 
particular instant. This can be done separately for each 
LCA eigenmode. As an analogon in TDLDA, we record 
the time evolution of the density n(r, t) and evaluate the 
Fourier components 



n(r,(jj^)= / n{r,t) eyip{—iuj^t)dt 



(6) 



for the frequencies lu^ that are associated with particu- 
lar excitations. Since this procedure is numerically more 
demanding, we have restricted the TDLDA analysis to 
a one dimensional section along the axis of symmetry, 
integrating over x and y coordinates. 

The top left picture in Fig. ^ shows a contour plot 
of the ground-state valence electron density of Na2 in 
a plane containing the axis of symmetry. (The grid for 
the calculation was larger than the shown part.) The 
ionic cores are clearly visible since they repel the va- 
lence electrons, leading to "holes" in the density. The 
bottom left picture visualizes how this valence electron 
density changes in time at the first excitation. Dark col- 
ors indicate a density increase, light colors a decrease. 
Obviously, the electron density increases at one end of 
the molecule and decreases at the other end. Thus, the 
valence electrons are shifted predominantly along the 
axis of symmetry. But the shift is not a uniform, sim- 
ple translation of the density along the dipole field (as it 
would be obtained from the sumrule estimate p6| ), but 
the intrinsic structure of the cluster is impressed on and 
reflected in the currents, leading to a shift "around" the 
ionic cores. The second excitation in LCA is a (twofold 
degenerate) x/y mode. Its density change (not shown 
in Fig. 1^) is predominantly perpendicular to the axis of 
symmetry, as one also naively would expect. The third 
LCA excitation, shown in the top right picture, is again 
a z mode. The regions of strongest density variation are 
shifted compared to the first z mode, and the oscilla- 
tion pattern shows a node at greater separation from the 
ionic cores. This reflects the mathematical requirement 
of orthogonality for the different modes fl|l7|,p7|. 
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Fig. 2 "Snapsliots" of the changes in valence electron den- 
sity associated with particular excitations of Na2. Contour 
plots show a plane containing the axis of symmetry of the 
molecule. Unit for the axes is the numerical grid spacing, 0.8 
ao- Top left: ground-state valence electron density n. Bot- 
tom left: dn/dt associated with the first excitation, i.e. first 
z-mode, in LCA. Top right: dn/dt associated with the third 
excitation, i.e. second z-mode, in LCA. Shadings lighter than 
the background gray indicate a density decrease, darker shad- 
ings an increase. Bottom right: n(r, lj^) integrated over x and 
y as a function of z for lowest (i.e., = 1) TDLDA excita- 
tion. The pictures are in accordance with understanding the 
excitations as density oscillations (see text). 



Physically, the plots show that the observed electronic 
transitions can be interpreted as different eigenmodes, 
i.e., intrinsic oscillation patterns of the valence electron 
distribution. The TDLDA transition densities confirm 
this picture. The bottom right part of Fig. |^ shows Eq. 
(^ evaluated for the frequency of the lowest mode. Since 
X and y coordinate have been integrated over, some of 
the finer structure might have been smoothed out. But 
it is clearly visible that also in TDLDA, the lowest exci- 
tation is associated with a density increase at one end of 
the molecule and a decrease at the other. The regions of 
maximum density change are found very similar in LCA 
and TDLDA. We have verified this also for the second 
excitation that is not shown in Fig. ^. Thus, the simple 
picture of excitations as density oscillations seems to be 
remarkably close to the truth, even for the small Na2. 

Fig. P shows the experimental low-temperature pho- 
toabsorption spectrum Q of Na^ and below the spec- 
tra obtained in TDSIC, TDLDA, and LCA. Again, the 
LCA results (bottom) are shown with a phenomenolog- 
ical line broadening to make comparison with TDLDA 
easier. Overall, the spectrum obtained in LCA is rather 
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Fig. 3 From top to bottom: Experimental [|], TDSIC, 
TDLDA and LCA photoabsorption spectrum S of Na^ in 
arbitrary units against excitation energy in eV. 



close to the experiment, which again is remarkable in 
view of the small cluster size. However, our main focus 
here is on the comparison between TDLDA and TDSIC. 
TDLDA gives the energies and relative peak heights for 
the two lower transitions close to the experimental ones. 
But instead of one peak that is seen experimentally at 
about 3.3 eV, TDLDA leads to two peaks at 3.20 eV 
and 3.53 eV. A similar pattern was also found but not 
explained in recent TDLDA calculations that fo- 
cused on the explanation of the observed linewidths. By 
going over to TDSIC, we can further investigate the na- 
ture of the double peak. The TDSIC spectrum in Fig. 
^ shows noticeable differences to TDLDA. First, a small 
subpeak is found at 0.92 eV. It is a relic of a Iph exci- 



tation which originally was close to 1 eV and which has 
given most of its strength to the dominant peak. The 
comparable state in TDLDA is found at about 1.5 eV, 
i.e., so close to the main peak that its strength is hardly 
recognizable. Second, whereas the peak at 2.09 eV is 
hardly changed by the averaged SIC, the peak which 
in TDLDA was at 2.76 eV shifts to 2.86 eV in TDSIC 
and appears broader since there is another transition 
close by at 3.03 eV. This might contribute to explaining 
why also in the experimental low temperature data, the 
middle peak appears to be somewhat broader. Finally, 
the last excitation again stays nearly unchanged at 3.57 
eV. Thus, by shifting the peak which in our TDLDA 
is found at 3.20 eV to lower energies, TDSIC leads to 
a spectrum that is close to the experimental one. We 
also tested whether this is only an indirect effect, due 
to slight rearrangements when the ionic geometry is re- 
optimized on SIC level. However, this is not the case: 
even when compared for exactly the same ionic struc- 
ture, TDLDA and TDSIC spectra show noticeable dif- 
ferences. We find, in accordance with earlier investiga- 
tions iQ, that TDSIC leaves the main resonance peaks 
basically unchanged, but it strongly modifies the single- 
particle energies, and thus the underlying Iph excita- 
tion spectrum. Our comparison with experimental data 
shows that while TDLDA gives reasonable results for 
the gross features of a spectrum, it can be inaccurate for 
details. In the energy range and for the clusters stud- 
ied here, the TDSIC description improves on TDLDA 
deficiencies in details of the coupling to Iph structures. 



4 Conclusions 

Our investigation of the photoabsorption spectra of two 
small sodium clusters with three different methods shed 
new light on the theoretical methods as well as on the 
understanding of the experiments. The local current ap- 
proximation is based on a "collective" picture of exci- 
tations. It exploits information that is contained in the 
curvature of the ground-state energy functional, and its 
success in the cases studied here demonstrates that the 
functional, indeed, contains relevant information on the 
excited states. Furthermore, LCA's quantitatively accu- 
rate description of the strong excitations and partial suc- 
cess in describing higher lying ones shows that the con- 
cept of "collectivity" and the detailed view in terms of 
particle-hole excitations have more in common than ex- 
pected, even for these small clusters. TDLDA has the 
advantage of being robust throughout a wide range of 
energies. It leads to a reliable description of overall fea- 
tures of photoabsorption spectra. The comparison with 
TDSIC showed, however, that details of the excitation 
spectrum can be rather sensitive to self-interaction ef- 
fects, even for a simple metal like sodium. Thus, treating 
exchange and correlation on a level beyond LDA is very 
desirable. 



6 



S. Kiimmel et al. 



The success of density-functional methods to accu- 
rately describe the excitation spectra of small clusters 
in general, and the LCA in particular, leads to an alter- 
native interpretation of the photoabsorption data. Our 
results show that the traditional way of thinking of ex- 
citations in small metal clusters as transitions between 
distinct molecular states coincides nicely with the more 
intuitive way of understanding them as collective elec- 
tronic eigenmodes, i.e., oscillations of the valence elec- 
tron density. Of course, these oscillations are neither ex- 
actly like the Mie plasmon in a classical metal sphere 
nor like the compressional bulk plasmon. Their frequen- 
cies and oscillation patterns are determined by the clus- 
ters' intrinsic structure, which for small systems like the 
ones studied here must of course be described quantum 
mechanically. But if this is taken into account, the pic- 
ture of density oscillations is well compatible with the 
"molecular states" point of view, and with experimental 
data. 
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